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Abstract 

The paper develops Newton's method of finding multiple eigenvalues with one 
Jordan block and corresponding generalized eigenvectors for matrices dependent on 
parameters. It computes the nearest value of a parameter vector with a matrix 
having a multiple eigenvalue of given multiplicity. The method also works in the 
whole matrix space (in the absence of parameters). The approach is based on the 
versal deformation theory for matrices. Numerical examples are given. 
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1 Introduction 



Transformation of a square nonsymmetric (non-Hermitian) matrix A to the Jordan canon- 
ical form is the classical subject that finds various applications in pure and applied mathe- 
matics and natural sciences. It is well known that a generic matrix has only simple eigen- 
values and its Jordan canonical form is a diagonal matrix. Nevertheless, multiple eigenval- 
ues typically appear in matrix families, and one Jordan block is the most typical Jordan 
structure of a multiple eigenvalue [31 E] . Many interesting and important phenomena as- 
sociated with qualitative changes in the dynamics of mechanical systems [2H1 EH1 EH1 EE] , 
stability optimization [31211 EH] and bifurcations of eigenvalues under matrix perturba- 
tions [33 ESI EH EH] are related to multiple eigenvalues. Recently, multiple eigenvalues with 
one Jordan block became of great interest in physics, including quantum mechanics and 
nuclear physics [21 El EI], optics [Sj, and electrical engineering [Hj. In most applications, 
multiple eigenvalues appear through the introduction of parameters. 

In the presence of multiple eigenvalues, the numerical problem of computation of the 
Jordan canonical form is unstable, since the degenerate structure can be destroyed by 
arbitrarily small perturbations (caused, for example, by round-off errors). Hence, instead 
of analyzing a single matrix, we should consider this problem in some neighborhood in 
matrix or parameter space. Such formulation leads to the important problem left open 
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by Wilkinson [JOl IS] : to find the distance of a given matrix to the nearest degenerate 
matrix. 

We study the problem of finding multiple eigenvalues for matrices dependent on several 
parameters. This implies that matrix perturbations are restricted to a specific submanifold 
in matrix space. Such restriction is the main difficulty and difference of this problem 
from the classical analysis in matrix spaces. Existing approaches for finding matrices 
with multiple eigenvalues [Zl El HH [T21 [TBI HH1 HH1 1231 12HI E31 SOI SU assume arbitrary 
perturbations of a matrix and, hence, they do not work for multiparameter matrix families. 
We also mention the topological method for the localization of double eigenvalues in two- 
parameter matrix families 

In this paper, we develop Newton's method for finding multiple eigenvalues with one 
Jordan block and corresponding generalized eigenvectors in multiparameter matrix fam- 
ilies. The presented method solves numerically the Wilkinson problem of finding the 
nearest matrix with a multiple eigenvalue (both in multiparameter and matrix space for- 
mulations). The implementation of the method in MATLAB code is available, see [3*T] . 
The method is based on the versal deformation theory for matrices. In spirit, our approach 
is similar to ^3], where matrices with multiple eigenvalues where found by path- following 
in matrix space (multiparameter case was not considered). 

The paper is organized as follows. In Section 2, we introduce concepts of singularity 
theory and describe a general idea of the paper. Section 3 provides expressions for values 
and derivatives of versal deformation functions, which are used in Newton's method in 
Section 4. Section 5 contains examples. In Section 6 we discuss convergence and accuracy 
of the method. Section 7 analyzes the relation of multiple eigenvalues with sensitivities of 
simple eigenvalues of perturbed matrices. In Conclusion, we summarize this contribution 
and discuss possible extensions of the method. Proofs are collected in the Appendix. 

2 Multiple eigenvalues with one Jordan block in mul- 
tiparameter matrix families 

Let us consider an m x m complex non-Hermitian matrix A, which is an analytical func- 
tion of a vector of complex parameters p = (pi, . . . ,p n ). Similarly, one can consider real 
or complex matrices smoothly dependent on real parameters, and we will comment the 
difference among these cases where appropriate. Our goal is to find the values of param- 
eter vector p at which the matrix A(p) has an eigenvalue A of algebraic multiplicity d 
with one d x d Jordan block (geometric multiplicity 1). Such and eigenvalue A is called 
nonderogatory. There is a Jordan chain of generalized vectors Ui , . . . , (the eigenvector 
and associated vectors) corresponding to A and determined by the equations 



Aui 
Au 2 



Aui, 

Au 2 + ui, 



(2.1) 



Au d 



Xu d + u d _i. 
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Figure 1: Geometry of the bifurcation diagram. 



These vectors form an m x d matrix U = [u 1; . . . , uj satisfying the equation 



AU = UJ, 



\ 
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A 
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1 

A/ 



(2.2) 



where 3\ is the Jordan block of size d. Recall that the Jordan chains taken for all the 
eigenvalues and Jordan blocks determine the transformation of the matrix A to the Jordan 
canonical form [T4] . 

In singularity theory parameter space is divided into a set of strata (smooth sub- 
manifolds of different dimensions), which correspond to different Jordan structures of the 
matrix A. Consider, for example the matrix family 



A(p) 
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P = {Pl,P2,Ps)- 



(2.3) 



The bifurcation diagram in parameter space is shown in Figure Q (for simplicity, we con- 
sider only real values of parameters). There are four degenerate strata: A 2 (surfaces), 
A 3 and Af A| (curves), and A 4 (a point). The surface A 2 , curve A 3 , and point A 4 corre- 
spond, respectively, to the matrices with double, triple, and quadruple eigenvalues with 
one Jordan block. The curve AfA 2 , is the transversal self-intersection of the stratum A 2 
corresponding to the matrices having two different double eigenvalues. This bifurcation 
diagram represents the well-known "swallow tail" singularity jl] . 

We study the set of parameter vectors, denoted by \ d , corresponding to matrices having 
multiple eigenvalues with one Jordan block of size d. The set X d is a smooth surface in 
parameter space having codimension d — 1 jSl E] • Thus, the problem of finding multiple 
eigenvalues in a matrix family is equivalent to finding the surface X d or its particular point. 
Since the surface X d is smooth, we can find it numerically by using Newton's method. This 
requires describing the surface \ d as a solution of d — 1 equations 



?i(p) = o, 



d. 



(2.4) 



3 



for independent smooth functions ?i(p). (In these notations, we keep the first function for 
the multiple eigenvalue A = <?i(p).) Finding the functions qi(p) and their first derivatives 
is the clue to the problem solution. 

In this paper, we define the functions qi(p) in the following way. According to versal 
deformation theory [311], in the neighborhood of X d , the matrix A(p) satisfies the relation 



/ ft(p) 



A(p)U(p) = U(p)B(p), B(p) 



\ 



92 IP) qi{P) 



(2.5) 



?i (P) / 
qd(p) are analytic func- 
, qd(p) are uniquely 



\ Qd(p) 

where U(p) is an m x d analytic matrix family, and gi(p), 
tions (blank places in the matrix are zeros). The functions qi(p), 
determined by the matrix family A(p). 

By using ()2.5jh it is straightforward to see that the surface X d is defined by equations 
(12. 4|) . If (|2.4|) are satisfied, the matrix B(p) is the dxd Jordan block. Hence, at p G X d , the 
multiple eigenvalue is A = qi(p) and the columns of U(p) are the generalized eigenvectors 
satisfying equations (|2.1jl . The method of finding the functions g*(p) and U(p) and their 
derivatives at the point p e X d has been developed in |2Tl I28j . In Newton's method for 
solving ()2.4|) . we need the values and derivatives of the functions qi(p) at an arbitrary 
point p ^ X d . 



3 Linearization of versal deformation functions 

Let po be a given parameter vector determining a matrix Ao = A(po). Since multi- 
ple eigenvalues are nongeneric, we typically deal with a diagonalizable matrix Ao. Let 
Ai, . . . , X m be eigenvalues of the matrix A . We sort these eigenvalues so that the first 
d of them, \i,...,\d, coalesce as the parameter vector is transferred continuously to the 
surface X d . The eigenvalues that form a multiple eigenvalue are usually known from the 
context of a particular problem. Otherwise, one can test different sets of d eigenvalues. 
Let us choose m x d matrices X and Y such that 

A X = XS, Y*A = SY*, Y*X = I, (3.1) 

where S is the dxd matrix whose eigenvalues are X\, . . . ,Xd] the star denotes the complex 
conjugate transpose. The first two equalities in (j3.1j) imply that the columns of the matrix 
X span the right invariant subspace of A corresponding to Ai, . . . , Xd, and the columns 
of Y span the left invariant subspace. The third equality is the normalization condition. 
The matrix S can be expressed as 

S = Y*A X, (3.2) 

which means that S is the restriction of the matrix operator Ao to the invariant subspace 
given by the columns of X. The constructive way of choosing the matrices S, X, and Y 
will be described in the next section. 
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The following theorem provides the values and derivatives of the functions ?i(p), . . . , <7d(p) 
in the versal deformation (|2.5|) at the point p . 

Theorem 3.1 Let S, Y , and X be the matrices satisfying equations \3.1\) . Then 

<7i(Po) = trace S/d, (3.3) 

and the values o/f^Po); • • • > Qd{Po) ore found as the characteristic polynomial coefficients 
of the traceless matrix S — gi(p )I.' 

z d - g 2 (po)^~ 2 gd-i(po)* - 9d(Po) = det ({z + gi(p ))I - S), (3.4) 

where I is the d x d identity matrix. The first derivatives of the functions qi(p) at po are 
determined by the recurrent formulae 



trace (y*?±7l) /d, 

((S - g 1 (p )I) J - 1 Y*^X^) - traceCC^ 1 )^- - ^trace^E 



klj 



i = 2, ...,d, j = l,...,n, 

(3-5) 

where the derivatives are evaluated at po/ C = B(po) — gi(po)I is the companion matrix 



Jq + X)gi(Po)Eii, (3-6) 



j=2 



/ 1 \ 

^(Po) ' ■ ■ 

V <?d(Po) o / 

and Eji z's t/ie matrix having the unit (i, l)th element and zeros in other places. 
The proof of this theorem is given in the Appendix. 

When the matrix A is arbitrary (not restricted to a multiparameter matrix family), 
each entry of the matrix can be considered is an independent parameter. Hence, the 
matrix A can be used instead of the parameter vector: p — > A. The derivative of A 
with respect to its (z, j)th entry is Ey. Thus, the formulae of Theorem 13. II can be applied. 

Corollary 3.1 Let S, Y, and X be the matrices satisfying equations \3.1\) . Then the 
values of qi(A ), . . . ,qd(A ) are given by formulae \3. 3]) and \3.$ with po substituted by 
A . Derivatives of the functions <7i(A), . . . , (^(A) with respect to components of the matrix 
A taken at A are 

|| = (XY7<i) T , 

|f = (X(S - aWfTf - trace(C")|f - £ trace^E^, (3 ' 7) 



Here T is the transpose operator, and 



( dqi 



dqj 
OA 



da 
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\ da m i 



dqi \ 

daim 

dqj 
da mm J 



(3i 



is the m x m matrix of derivatives of qi(A) with respect to components of the matrix A 
taken at A . 

At p G X d , we can find the multiple eigenvalue A and the corresponding Jordan chain 
of generalized eigenvectors u 1; . . . , u^. This problem reduces to the transformation of the 
matrix S to the prescribed Jordan form (one Jordan block). A possible way of solving this 
problem is presented in the following theorem (see the Appendix for the proof.). 

Theorem 3.2 At the point p G \ d , the multiple eigenvalue is given by the expression 

A = traceS/d. (3.9) 

The general form of the Jordan chain of generalized eigenvectors ui, . . . , is 

X(S-AI) d - 1 k, u d _! = X(S - AI)k, u d = Xk, (3.10) 



Ui 



where k G C d is an arbitrary vector such that the eigenvector ui is nonzero. Choosing a 
particular unit-norm eigenvector u±, e.g., by taking the scaled biggest norm column of the 
matrix X(S — AI) d_1 ; one can fix the vector k by the orthonormality conditions 



d. 



(3-11) 



The accuracy of the multiple eigenvalue and generalized eigenvectors determined by 
formulae \3. ty) and has the same order as the accuracy of the point po in the 

surface \ d . 



4 Newton's method 

There are several ways to find the matrices S, X, and Y. The simplest way is to use the 
diagonalization of Ao- Then S = diag(Ai, . . . , A^) is the diagonal matrix, and the columns 
of X and Y are the right and left eigenvectors corresponding to Ai, . . . , A^. This way will 
be discussed in Section 5. 

If the parameter vector po is close to the surface \ d , the diagonalization of the matrix 
Ao is ill-conditioned. Instead of the diagonalization, one can use the numerically sta- 
ble Schur decomposition S = X*AoX, where S is an upper-triangular matrix called the 
Schur canonical form, and X = (X*)" 1 is a unitary matrix [15J. The diagonal elements 



6 



Sii, . . . , s mm of S are the eigenvalues of Ao- We can choose the Schur form so that the first 
d diagonal elements are the eigenvalues Ai, . . . , A^. Performing the block-diagonalization 
of the Schur form S [To"! §7.6], we obtain the block-diagonal matrix 

(o S') = [Y,Y']*A [X,X], (4.1) 

where S is a d x d upper-triangular matrix with the diagonal (Ai, . . . , Ad); [X, X'] and 
[Y, Y']* = [X, X'] -1 are nonsingular m x m matrices (not necessarily unitary). These 
operations with a Schur canonical form are standard and included in many numerical 
linear algebra packages, for example, LAPACK pQ. They are numerically stable if the 
eigenvalues Ai, . . . , A^ are separated from the remaining part of the spectrum. As a result, 
we obtain the matrices S, X, and Y satisfying equations QH.lj) . 

When the matrices S, X, and Y are determined, Theorem 13. II provides the necessary 
information for using Newton's method for determining the stratum X d . Indeed, having 
the parameter vector p = (p°, • • • ,P°) as the initial guess, we linearize equations (|2.4j) of 
the surface X d as 

n f> 

QiM+^2j L {Pi-Pj)=^ i = 2,-..,d, (4.2) 

3=1 3 

where the values of ft(po) and the derivatives dqi/dpj at po are provided by Theorem 13.11 
In the generic case, the linear part in (|4.2|) is given by the maximal rank matrix [dqi/dpj]. 
System (|4.2|) has the single solution if the number of parameters n = d — 1 (the set X d 
is an isolated point). If n > d — 1, one can take the least squares solution or any other 
solution depending on which point of the surface X d one would like to find. If n < d, 
the multiple eigenvalue still can exist in matrices with symmetries (e.g., Hamiltonian or 
reversible matrices [35]); then the least squares fit solution of (|4.2|) is a good choice. 

In Newton's method, the obtained vector of parameters p = (pi, . . . ,p n ) is used in the 
next iteration. In each iteration, we should choose d eigenvalues of the matrix A . These 
are the d eigenvalues nearest to the approximate multiple eigenvalue 

n F) 

A = <&(pd) + £^(p;-pS) (4 ' 3) 

j=i 

calculated at the previous step. If the iteration procedure converges, we obtain a point 
p G X d . Then the multiple eigenvalue and corresponding generalized eigenvectors are 
found by Theorem 13.21 Note that, at the point p G \ d , system (|4.2|) determines the 
tangent plane to the surface X d in parameter space. The pseudo-code of the described 
iteration procedure is presented in Table ^ Depending on a particular application, the 
line 3 in this pseudo-code can be implemented in different ways, e.g., as the least squares 
solution or as the solution nearest to the input parameter vector p . The implementation 
of this method in MATLAB code is available, see |31j . 

In case of complex matrices dependent on real parameters, the same formulae can be 
used. In this case, system (14.2(1 represents 2{d — 1) independent equations (each equality 
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INPUT: matrix family A(p), initial parameter vector po, and eigenvalues Ai, . . . , Ad 



Schur decomposition and block-diagonalization f)4. 1|) of the matrix Ao = A(po); 

evaluate %(po) and dqi/dpj by formulae (|3.3J1 - (|3.5|) : 

find Pnew by solving system (|4.2jl (e.g. the least squares solution); 

IF 1 1 p„e«; — Poll > desired accuracy 

evaluate approximate multiple eigenvalue \ app by (j4.3jl : 

choose d eigenvalues A™ 6 ™, . . . , \% ew of A new = A(p new ) nearest to A app ; 

perform a new iteration with p = p n ew and Aj = A™ 6 "", i — 1, . . . , d (GOTO 1); 
ELSE (IF \\pnew — Poll < desired accuracy) 

find multiple eigenvalue and generalized eigenvectors by formulae (|3.9|) - (|3.11|) ; 



OUTPUT: parameter vector p £ A d , multiple eigenvalue A and Jordan chain 
of generalized eigenvectors u 1; . . . , u d 

Table 1: Pseudo-code of Newton's method for finding multiple eigenvalues in multiparam- 
eter matrix families. 

determines two equations for real and imaginary parts). This agrees with the fact that 
the codimension of X d in the space of real parameters is 2(d — 1) 

Finally, consider real matrices smoothly dependent on real parameters. For complex 
multiple eigenvalues, the system (|4.2jl contains 2(d — 1) independent real equations (codi- 
mension of X d is 2(d - 1)). Remark that imaginary parts of the eigenvalues Ai, . . . , A^ 
should have the same sign. For real multiple eigenvalues, (&(po) and dqi/dpj are real 
(the real Schur decomposition must be used). Hence, ()4.2j) contains d — 1 real equations 
(codimension of X d is d — 1). In this case, the eigenvalues Ai, . . . , A^ are real or appear in 
complex conjugate pairs. 

In some applications, like stability theory [35j, we are interested in specific multiple 
eigenvalues, e.g., zero and purely imaginary eigenvalues. In this case equation (|4.3|) should 
be included in the linear system of Newton's approximation ([4.2)1 . 

For arbitrary matrices A = [ay] (without parameters), similar Newton's iteration 
procedure is based on Corollary 13.11 The linearized equations f)4.2p are substituted by 



dqi 

j,k=i 



&(po) + d^i^ ~ a ° fc) = °' % = 2 ' • • • ' d ' (4 - 4) 



where A = [a° fe ] is the matrix obtained at the previous step or the initial input matrix. 
The first-order approximation of the multiple eigenvalue (|4.3J) takes the form 

A = ?i(po) + d^~^ ajk ~ a °^- ^ 

j,k=l Jk 
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5 Examples 



All calculations in the following examples were performed using MATLAB code [31 . For 
the sake of brevity, we will show only first several digits of the computation results. 

5.1 Example 1 

Let us consider the two-parameter family of real matrices 



A(p) = p 1 1 p 2 , p=(pi,p 2 ). (5.1) 




Bifurcation diagram for this matrix family is found analytically by studying the discrimi- 
nant of the characteristic polynomial. There are two smooth curves A 2 and a point A 3 at 
the origin (the cusp singularity), see Figure El 

Let us consider the point p = (—0.03, 8.99), where the matrix A has the eigenvalues 
Ai 2 = —1.995 ± «0.183 and A 3 = 6.990. In order to detect a double real eigenvalue, we 
choose the pair of complex conjugate eigenvalues Ai, A 2 . By ordering diagonal blocks in 
the real Schur form of A and block-diagonalizing, we find the matrices S, X, and Y 
satisfying (|3.1|) in the form 



-1.995 -5.083 
0.007 -1.995 



0.688 -0.676 \ / 0.729 -0.574 

X = I -0.688 -0.491 I , Y = | -0.604 -0.286 



0.231 0.550 / \ 0.357 0.858 

Applying the formulae of Theorem 13.11 we find 

fc(po) \ = ( -1-995 \ / d qi /d Pl dqi/dp 2 \ = / -0.111 -0.148 
g 2 (po) J V -°- 033 J ' V dq 2 /dp 1 dq 2 /dp 2 J V L001 °- 333 



(5.2) 



The linearized system (14. 2|) represents one real scalar equation. We find the nearest pa- 
rameter vector p G A 2 (the least squares solution) as 

P = Po - jr. m m ^2 ( 9( l2/dpi, dq 2 /dp 2 ) = (-0.00001, 8.99999). (5.3) 

{dq 2 /dpiY + (dq 2 /dp 2 y 

After five iterations of Newton's method, we find the exact nearest point p = (0, 9) G A 2 . 
Then Theorem 13.21 gives the multiple eigenvalue and the Jordan chain with the accuracy 
10- 15 : 

j / 3 -1 + 30/19 \ 
A = -2, [ Ul ,u 2 ] = ^ -3 2 -30/19 . (5.4) 
^ 19 \ 1 -1 + 10/19 / 

Now let us take different points po in the neighborhood of the curve A 2 and calculate 
one-step Newton's approximations of the nearest points p G A 2 . In this case we choose 
Ai, A 2 as a pair of complex conjugate eigenvalues of A = A(p ). If all eigenvalues of A 
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Figure 2: One-step approximations of the nearest points with double eigenvalues. 



are real, we test all different pairs of eigenvalues, and take the pair providing the nearest 
point p G A 2 . The result is shown in Figure El where each arrow connects the initial 
point po with the one-step Newton's approximation p. For one point po we performed 
two iterations, taking the point p as a new initial point po = p. The convergence of 
this iteration series is shown in the enlarged part of parameter space (inside the circle in 
Figure |2j). The results confirm Newton's method rate of convergence. 



5.2 Example 2 

Let us consider the real matrix A = Ai + eE, where 





i 
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A x = 8 , E = 8 3 6 , (5.5) 



and £ = 2.2e-15, 5 = 1.5e-9. This matrix was used in HO] for testing the GUPTRI PjlTHj 
algorithm. It turned out that this algorithm detects a matrix A G A 3 (with a nonderoga- 
tory triple eigenvalue) at the distance O(10 -6 ) from Ao, while the distance from Ao to 
A 3 is less than ||eE|| F = 3.62e-14 since Ai G A 3 . This is explained by the observation 
that the algorithm finds matrix perturbations along a specific set of directions, and these 
directions are almost tangent to the stratum A 3 in the case under consideration [TT]] . 

Our method determines locally the whole stratum A 3 in matrix space and, hence, it 
should work correctly in this case. Since the triple eigenvalue is formed by all eigenvalues 
of A Q , we can use S = A and X = Y = I in the formulae of Corollary 13.11 As a result, 
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we find the least squares solution of system ()4.4j) in the form A = Aq + AA, where 



/ \ 
AA=1.0e-14* -1.760 0. (5.6) 
\ -0.880 / 

Approximations of the multiple eigenvalue and corresponding generalized eigenvectors 
evaluated by Theorem 13. 21 for the matrix A are 

/ 1.000 -0.000 -0.000 \ 
A = 8.800e-15, [u 1; u 2 , u 3 ] = 0.000 1.000 -0.000 . (5.7) 

\ 0.000 0.000 6.667e+8 / 

We detected the matrix A at the distance ||AA||^ = 1.97e-14, which is smaller than 
the initial perturbation ||£E||^ = 3.62e-14 (||AA||^ denotes the Frobenius matrix norm). 
The matrix U = [111,112,113] satisfies the Jordan chain equation ()2.2|) with the very high 
accuracy || AU - UJ a ||f/||U|| f = 9.6e-23. 

The normal complementary subspace N of the tangent space to A 3 at A x has the 
form [3] 

f /o 0\ , 

N = < \ x I x,y e M>. (5.8) 
L \ y 5x / > 

It is easy to see that the matrix AA in (J5.6|) is equal to the projection of — eE to the 
normal subspace N. This confirms that the obtained matrix A e A 3 is the nearest to A . 

5.3 Example 3 

Let us consider the 12 x 12 Frank matrix Ao = [a° ] with the elements 

„o / n + l-max(i,j), j>i-l, ft . Q s 
«y -| , j < i — 1. (5 ' 9) 

The Frank matrix has six small positive eigenvalues which are ill-conditioned and form 
nonderogatory multiple eigenvalues of multiplicities d = 2, . . . , 6 for small perturbations 
of the matrix. The results obtained by Newton's method with the use of Corollary 13.11 
are presented in Table 121 An eigenvalue of multiplicity d of the nearest matrix A e X d is 
formed by d smallest eigenvalues of Ao- The second column of Table 121 gives the distance 
dist(A , X d ) = || A — Ao || f, where the matrix A is computed after one step of Newton's 
procedure. The third column provides exact distances computed by Newton's method, 
which requires 4-5 iterations to find the distance with the accuracy O(10 -15 ). At each 
iteration, we find the solution A of system (14 .4j) . which is the nearest to the matrix 
(15. 9|) . The multiple eigenvalues and corresponding generalized eigenvectors are found at 
the last iteration by Theorem 13 .21 The accuracy estimated as || AU — UJa||f/||U||f varies 
between 10 -10 and 10~ 13 . The matrices of generalized eigenvectors U have small condition 
numbers, which are given in the fourth column of Table 121 For comparison, the fifth and 
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d 


dist(A ,A d ) 
1-step approximation 


dist(A ,A d ) 
exact 


condU 


||A- A ||f Q3] 


||A- A ||f ca 


2 


1.619e-10 


1.850e-10 


1.125 


3.682e-10 




3 


1.956e-8 


2.267e-8 


1.746 


3.833e-8 




4 


1.647e-6 


1.861e-6 


4.353 


3.900e-6 




5 


9.299e-5 


1.020e-4 


14.14 


4.280e-4 


6e-3 


6 


3.150e-3 


3.400e-3 


56.02 


7.338e-2 





Table 2: Distances to the multiple eigenvalue strata A for the Frank matrix, 
sixth columns give upper bounds for the distance to the nearest matrix A £ X d found 

in pang. 

We emphasize that this is the first numerical method that is able to find exact distance 
to a nonderogatory stratum X d . Methods available in the literature cannot solve this 
problem neither in matrix space nor for multiparameter matrix families. 

6 Convergence and accuracy 

In the proposed approach, the standard Schur decomposition and block-diagonalization 
(14. lj) of a matrix are required at each iteration step. Additionally, first derivatives of the 
matrix with respect to parameters are needed at each step. Numerical accuracy of the 
block-diagonalization depends on the separation sep(S, S') of the diagonal blocks in the 
Schur canonical form (calculated prior the block-diagonalization) Instability occurs 
for very small values of sep(S,S'), which indicates that the spectra of S and S' overlap 
under a very small perturbation of A . Thus, numerical instability signals that the chosen 
set of d eigenvalues should be changed such that the matrix S includes all "interacting" 
eigenvalues. 

The functions qi(p) are strongly nonlinear near the boundary of the surface X d . The 
boundary corresponds to higher codimension strata associated with eigenvalues of higher 
multiplicity (or eigenvalues of the same multiplicity but with several Jordan blocks). For 
example, the stratum A 2 in FigureQis bounded by the singularities A 3 and A 4 . As a result, 
the convergence of Newton's method may be poor near the boundary of X d . This instability 
signals that we should look for eigenvalues with a more degenerate Jordan structure (e.g. 
higher multiplicity d) . Analysis of the surface X d very close to the boundary is still possible, 
but the higher precision arithmetics may be necessary. 

Figure |H] shows first iterations of Newton's procedure for different initial points in 
parameter space for matrix family (J5.1|) from Example 1. Solid arrows locate double 
eigenvalues (the stratum A 2 ) and the dashed arrows correspond to triple eigenvalues (the 
stratum A 3 ). One can see that the stratum A 2 is well approximated when po is relatively far 
from the singularity (from the more degenerate stratum A 3 ). For the left-most point po in 
Figure EH the nearest point p £ A 2 simply does not exist (infimum of the distance ||p — p || 
for p £ A 2 corresponds to the origin p = £ A 3 ). Note that, having the information on 
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Figure 3: One-step approximations of the nearest points of the strata A 2 (solid arrows) 
and A 3 (dashed arrows) near the cusp singularity. 

the stratum A 3 , it is possible to determine locally the bifurcation diagram (describe the 
geometry of the cusp singularity in parameter space) (23 EH] • 

For the backward error analysis of numerical eigenvalue problems based on the study 
of the pseudo-spectrum we refer to [37|. We note that the classical numerical eigenvalue 
problem is ill-conditioned in the presence of multiple eigenvalues. The reason for that is the 
nonsmoothness of eigenvalues at multiple points giving rise to singular perturbation terms 
of order £ 1//d , where d is the size of Jordan block [35]. On the contrary, in our problem we 
deal with the regular smooth objects: the strata X d and the versal deformation B(p). 

7 Approximations based on diagonal decomposition 

In this section we consider approximations derived by using the diagonal decomposition 
of Ao- The diagonal decomposition is known to be ill-conditioned for nearly defective 
matrices. However, this way is easy to implement, while the very high accuracy may be 
not necessary. According to bifurcation theory for eigenvalues |3S1, the accuracy of the 
results based on the diagonal decomposition will be of order e 1 ^, where e is the arithmetics 
precision. Another reason is theoretical. Bifurcation theory describes the collapse of 
a Jordan block into simple eigenvalues pffi EE]- Our approximations based on the 
diagonal decomposition solve the inverse problem: using simple (perturbed) eigenvalues 
and corresponding eigenvectors, we approximate the stratum X d at which these eigenvalues 
coalesce. 

Let us assume that the matrix Ao is diagonalizable (its eigenvalues Ai,...,A m are 
distinct). The right and left eigenvectors of A are determined by the equations 

A Xi = AjXi, y*A = Xiy*, y* x » = 1 I 7 - 1 ) 

with the last equality being the normalization condition. 

In Theorem 13.11 we take S = diag(Ai, . . . , A^) = Y*A X, X = [xi, . . . , xj, and 
Y = [yi,...,yj, where X 1 ,...,X d are the eigenvalues coalescing at p G X d . In this 
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case expressions ()3.5)1 take the form 
dqi = 1 y> „gA 



^ = X> - gi(Po))^ 1 y^x, - tracelC- 1 )^ - gtrace(C^E fel )^, 
d Pj d Pj d Pj f^ 2 d Pj 

i = 2, j = 1, . . . ,n. 

The interesting feature of these expressions is that they depend only on the simple eigen- 
values Ai, . . . , Xd and their derivatives with respect to parameters at po [35] : 

Wj =Yi W *> . = l,...,d. (7.3) 

For example, for d = 2 we obtain the first-order approximation of the surface A 2 in the 
form of one linear equation 

*+'E(*g*-y^)fe-i*> = o, * = ^- (7 - 4) 

Let us introduce the gradient vectors VAj = (dXi/dpi, . . . , d\i/dp n ), i = 1,2, with the 
derivatives dXi/dpj given by expression (|7.3|) . Then the solution of (J7.4|) approximating 
the vector p G A 2 nearest to p is found as 



VA2-VA1 

P = P0 " ||VA 2 -VA 1 ||^ - (7 ' 5) 

It is instructive to compare this result with the first-order approximation of the nearest 
p G A 2 , if we consider Ai and A2 as smooth functions 

A J ( P ) = A l + ^^fe-^) + 0(|| P - P oH 2 ), < = 1,2. (7.6) 

j=i 

Using ()7.6|) in the equation Ai(p) = A 2 (p) and neglecting higher order terms, we find 
which yields the nearest p as 



VA 2 - VAi nx 
P = Po ~ ||VA 2 - VA.P 2 ' 5 - (7 ' 8) 

Comparing (|7.8j) with (|7.5j) . we see that considering simple eigenvalues as smooth 
functions, we find the correct direction to the nearest point p G A 2 , but make a mistake 
in the distance to the stratum A 2 overestimating it exactly twice. This is the consequence 
of the bifurcation taking place at p G A 2 and resulting in 0(||p — Poll 1 ^ 2 ) perturbation of 
eigenvalues and eigenvectors [321 HHH EE! • 
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8 Conclusion 



In the paper, we developed Newton's method for finding multiple eigenvalues with one 
Jordan block in multiparameter matrix families. The method provides the nearest pa- 
rameter vector with a matrix possessing an eigenvalue of given multiplicity. It also gives 
the generalized eigenvectors and describes the local structure (tangent plane) of the stra- 
tum X d . The motivation of the problem comes from applications, where matrices describe 
behavior of a system depending on several parameters. 

The whole matrix space has been considered as a particular case, when all entries of a 
matrix are independent parameters. Then the method provides an algorithm for solving 
the Wilkinson problem of finding the distance to the nearest degenerate matrix. 

Only multiple eigenvalues with one Jordan block have been studied. Note that the 
versal deformation is not universal for multiple eigenvalues with several Jordan blocks 
(the functions <?i(p), . . . , <?d(p) are n °t uniquely determined by the matrix family) Hj. 
This requires modification of the method. Analysis of this case is the topic for further 
investigation. 

9 Appendix 

9.1 Proof of Theorem fXTI 

Taking equation (|2.5j) at po, we obtain 

A U = U B , (9.1) 

where U = U(p ) and B = B(p ). Comparing ()9.1|) with ([3.1)1 . we find that the matrix 
B is equivalent up to a change of basis to the matrix S. Then the equality ([3.3)1 is 
obtained by equating the traces of the matrices Bo and S, where Bo has the form ([2.5)1 . 
Similarly, the equality ([3.4jl is obtained by equating the characteristic equations of the 
matrices B — gi(po)I and S — gi(po)I. 

The columns of the matrices X and Uo span the same invariant subspace of A . Hence, 
the matrices X and Uo are related by the expression 

Uo = XF, (9.2) 

for some nonsingular d x d matrix F. Using ([3.1jl and ([9.2)1 in ([9.1)1 . we find the relation 

SF = FB . (9.3) 

Taking derivative of equation ([2.5)1 with respect to parameter pj at p , we obtain 

t <9U <9U„ OB <9A TT . , 

OPj Opj Opj Opj 

Let us multiply both sides of ([9.4)1 by the matrix F _1 (S — gi(po)I) 1-1 Y* and take the 
trace. Using expressions ([3.1)1 . ([9.3)1 . and the property trace(AB) = trace(BA), it is 
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straightforward to check that the left-hand side vanishes and we obtain the equation 

= trace ((u g - f^F^S - ^(poW^Y*) . (9.5) 

Substituting (|9.2j) into ()9.5|) and using equalities (j3.1|) . ()9.3|) . and B = g!(p )I + C, we 
find 

trace ( C^ 1 ^ ] - trace ( (S - gi(p )I) <_1 Y*^X ] = 0. (9.6) 
V dPjJ V d Pj J 

Using ()2.5|) . ()3.6|) and taking into account that trace(C* _1 Efci) = for k > i, we obtain 



trace(C i - 1 )^ + V trace(C^ 1 E fcl )^ = trace ( (S - 9l (p 



jNi-ly^X . (9.7) 

% / 

Taking equation (|9.7|) for i = 1, . . . , d and using the equality trace(C*~ 1 E il ) = 1, we get 
the recurrent procedure ()3.5|) for calculation of derivatives of (?i(p), . . . , <?d(p) at p . 

9.2 Proof of Theorem 15^1 

At po G X d , we have ^(po) = • • • = ^(po) — 0. From (|2.5j) it follows that B = Ja is 
the Jordan block with the eigenvalue A = gi(po) = trace S/d. By using (|3.1j) . one can 
check that the vectors ()3.1()j) satisfy the Jordan chain equations (j2.1)l for any vector k. 
Equations (j3.11j) provide the way of choosing a particular value of the vector k. 

Since B = Ja, the versal deformation equation ()2.5|) becomes the Jordan chain equa- 
tion ()2.2|) at po £ X d - Hence, the columns of the matrix Uo are the generalized eigenvectors. 
Since the function qi(p) and the matrix U(p) smoothly depend on parameters, the ac- 
curacy of the multiple eigenvalue and generalized eigenvectors has the same order as the 
accuracy of p G X d . 
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